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Two  different  types  of  simulations  are  required  for  analyzing  Army  vehicles  when  they  are  subjected  to  explosive/impact  loads:  (a) 
Survivability  analysis,  when  loads  generate  damage  on  the  vehicle  structure  by  causing  permanent  deformation;  (b)  Shock  analysis,  when 
the  vehicle  structure  remains  within  the  elastic  region.  Shock  analysis  ensures  that  vibration  levels  must  remain  low  at  locations  where 
electronic  equipment  is  mounted.  Due  to  short  duration  of  the  load  the  high  frequency  response  is  important  for  shock  analysis  (i.e.  500Hz 
-  10,000Hz).  Since  the  shock  loads  comprise  one  of  the  sets  of  design  loads  for  any  Army  vehicle,  it  is  important  to  perform  a  shock 
simulation  efficiently  enough  in  order  to  make  design  decisions.  Such  simulations  will  allow  to  access  the  survivability  of  the  equipment 
and  therefore  of  the  vehicle,  and  they  will  allow  to  incorporate  design  changes  in  order  to  achieve  desirable  shock  spectra  at  mounting 
locations.  Computing  efficiently  the  high  frequency  response  of  composite  structures  is  the  objective  of  this  research. 
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Modeling  High  Frequency  Vibration  in  Composites  Using  an  Energy  Finite  Element 
Method  for  Shock  Analysis  of  Lightweight  Army  Vehicles 

(Proposal  No.  51087EG) 

Nickolas  Vlahopoulos 
Professor 

Department  of  Naval  Architecture  and  Marine  Engineering 
Department  of  Mechanical  Engineering 
University  of  Michigan,  Ann  Arbor,  MI  48109 

I.  Research  Objectives 

•  Two  different  types  of  simulations  are  required  for  analyzing  Army  vehicles  when  they 
are  subjected  to  explosive/impact  loads. 

(a)  Survivability  analysis;  when  loads  generate  damage  on  the  vehicle  structure  by  causing 
permanent  deformation  (typical  results  are  presented  in  Figure  1). 


Figure  1.  Survivability  analysis  for  vehicle/occupant  response  to  explosive  threat  [N. 
Vlahopoulos,  G.  Zhang,  “Validation  of  a  Simulation  Process  for  Assessing  the  Response  of  a 
Vehicle  and  its  Occupants  to  an  Explosive  Threat,”  27th  Army  Science  Conference] 


(b)  Shock  analysis;  when  the  vehicle  structure  remains  within  the  elastic  region.  Shock  analysis 
ensures  that  vibration  levels  must  remain  low  at  locations  where  electronic  equipment  is 
mounted.  Due  to  short  duration  of  the  load  the  high  frequency  response  is  important  for  shock 
analysis  (i.e.  500Hz  -  10,000Hz).  Typical  shock  analysis  results  are  presented  in  Figure  2. 
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Figure  2.  (a)  Plate  assembly  structure  analyzed  for  shock;  (b)  response  spectra  at  a  point 
[R.  D.  Hampton,  N.  S.  Wiedenman,  T.  H.  Li,  “Analytical  Determination  of  Shock  Response 
Spectra,”  1 1th  ARL/USMA  Technical  Symposium] 

•  In  an  effort  to  make  Army  vehicles  more  lightweight,  composite  materials  can  be  used 
for  their  construction.  Efficient  shock  vibration  analysis  of  composite  vehicle  comprises 
the  objective  of  this  research  effort. 
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II.  Approach 

•  Traditionally,  Finite  Element  Analysis  (FEA)  solutions  are  used  in  order  to  compute  the 
vibration  of  a  structure  due  to  blast/shock.  These  conventional  FEA  methods  compute 
the  displacement  of  the  vibration  at  each  point  of  the  structure  with  respect  to  time.  At 
high  frequencies  the  wavelengths  of  the  vibration  are  small  and  a  large  number  of  nodes 
and  elements  are  required  in  a  FEA  model  in  order  to  capture  correctly  the  shape  of  the 
vibration  due  to  the  small  wavelengths.  The  large  number  of  elements  in  the  FEA  model 
result  in  requirements  for  extremely  high  computational  resources.  Thus,  high  frequency 
shock  analyses  for  an  entire  vehicle  using  conventional  FEA  methods  are  very  expensive 
and  even  infeasible  after  a  certain  frequency. 

•  The  Energy  Finite  Element  Analysis  (EFEA)  formulation  is  based  on  using  the  space 
averaged  over  a  wavelength  and  time  averaged  over  a  period  energy  density  (energy 
density)  as  a  primary  variable.  The  governing  differential  equations  of  the  EFEA 
formulation  are  developed  for  the  energy  density.  The  energy  density  varies 
exponentially  in  space,  thus  only  a  small  number  of  elements  is  required  in  order  to 
capture  the  high  frequency  vibration.  Figure  3  depicts  graphically  the  differences  in  the 
primary  variables  between  conventional  FEA  and  EFEA  methods.  Once  the  governing 
differential  equations  are  identified,  the  finite  element  approach  can  be  employed  for 
generating  the  corresponding  numerical  system  of  equations  and  the  element  matrices. 
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Figure  3.  Graphical  representation  of  the  differences  between  conventional  FEA  and 
EFEA  and  the  main  reason  for  efficiency  in  EFEA  calculations  due  to  the  small  number 

of  elements 


The  main  challenges  in  the  EFEA  arise  when  different  structural  components  are 
connected  together.  Unlike  the  conventional  FEA,  where  the  primary  variable  remains 
continuous  at  the  connection,  in  the  EFEA  the  primary  variable  becomes  discontinuous. 
This  situation  is  represented  graphically  in  Figure  4.  Specialized  joint  formulations  must 
be  developed  for  connecting  the  element  level  matrices  across  a  joint  based  on  power 
balance  equations. 

Energy  density  discontinuous  across  a  joint 

1 


Displacement  continuous  across  a  joint 

Figure  4.  Graphical  representation  of  the  challenges  in  the  EFEA  at  connections 

between  components 

In  order  to  develop  an  EFEA  formulation  for  composite  structures  two  main  new 
developments  must  be  completed  and  implemented  into  software: 

(i)  Development  of  the  governing  differential  equations  for  composite  materials  and  the 
associated  element  matrices. 

(ii)  Development  of  joint  formulations  for  the  connections  between  composite 
components. 

These  two  tasks  comprise  the  main  focus  of  this  research  project. 


Significance 


In  a  continuous  effort  to  reduce  fuel  consumption  of  Army  vehicles,  light  weight  vehicle 
structures  with  increased  functionality  must  be  designed.  Since  the  shock  loads  comprise 
one  of  the  sets  of  design  loads  for  any  Army  vehicle,  it  is  important  to  perform  a  shock 
simulation  efficiently  enough  in  order  to  make  design  decisions.  Such  simulations  will 
allow  to  access  the  survivability  of  the  equipment  and  therefore  of  the  vehicle,  and  they 
will  allow  to  incorporate  design  changes  in  order  to  achieve  desirable  shock  spectra  at 
mounting  locations. 

Currently,  typical  computational  times  based  on  conventional  finite  element  methods 
require  two  days  per  impact  load  (Figure  2).  From  the  response  spectra  it  can  be 
observed  that  the  majority  of  the  shock  energy  appears  in  the  frequency  range  1,000Hz  - 
10,000Hz. 

The  EFEA  is  an  established  simulation  method  for  metallic  structures  (Figure  5).  In  this 
project  research  has  been  conducted  for  enabling  the  utilization  of  the  EFEA  for 
analyzing  composite  structures. 


Figure  5.  EFEA  applications  for  metallic  structures 


Accomplishments  in  the  current  grant 

Development  of  the  governing  differential  equations  and  element  level  matrices  for 
composite  materials.  Appendix  A  contains  detailed  information  about  the  derivation  of 
the  element  level  EFEA  matrices  for  composite  materials. 

Development  of  joint  formulations  for  the  connections  between  composite  components. 
Appendix  B  contains  technical  information  about  the  computation  of  the  power  transfer 
coefficients  and  of  the  joint  matrices. 

Use  the  element  matrices  and  the  joint  matrices  for  creating  the  system  of  equations  that 
represents  complex  composite  structures 


•  Gradual  validation  through  comparison  with  conventional  finite  element  solutions  for 
simple  structures  (Figures  6  and  7)  and  for  a  complex  structure  (Figures  8  and  9)  similar 
to  the  one  tested  under  explosive  loads  in  Figure  1 .  For  the  composite  structure  with  a 
double  floor  and  outer  V  shaped  bottom  presented  in  Figure  8  the  excitation  is  applied  in 
the  middle  of  the  outer  V  shaped  floor. 


The  model  in  conventional  FE A  has  12.800  elements. 

The  model  in  EFE  A  has  only  32  elements. 

Figure  6.  Conventional  FEA  and  EFEA  models  for  two  composite  panels 
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Computation  time:  6.5  hours  vs.  5  minutes 

Figure  7.  Comparison  between  FEA  and  EFEA  results  for  two  different  configurations  of 
connecting  composite  panels  one  and  two  along  two  different  edges  of  panel  1 


Figure  8.  All  composite  structure  with  a  double  floor  and  outer  V  shaped  bottom;  EFEA  model 


Figure  9.  Comparison  between  conventional  FEA  and  EFEA  results  for  the  response  of  the  various  panels 

Figure  9  presents  the  correlation  between  the  conventional  FEA  model  (established  computational 
method)  and  the  new  EFEA  method  for  modeling  composite  structures.  The  agreement  between  the  two 
solutions  in  terms  of  capturing  the  frequency  dependency  of  the  vibration  at  the  different  parts  of  the 
structure,  along  with  the  magnitude  of  the  response  demonstrate  the  validity  of  using  the  EFEA  for 
modeling  the  transfer  of  vibrational  energy  in  composite  structures.  The  analysis  was  terminated  at 
1,000Hz  because  this  is  the  upper  limit  of  validity  of  the  conventional  FEA  model  used  in  this  study. 
There  is  no  upper  frequency  limit  in  the  EFEA  model.  In  this  case  study  the  EFEA  simulations  are  25 
times  faster  compared  to  the  FEA  for  the  same  desktop  computer. 

V.  Technology  Transfer 

Continuous  interactions  is  in  place  with  Farzad  Rostam  -Abedi  and  Kris  Bishnoi  from  the  survivability 
group  of  TARDEC.  Visits  and  seminars  on  this  topic  have  been  given  at  TARDEC.  Several  interactions 
have  also  taken  place  with  Matt  Castanier  from  the  CASI  group  of  TARDEC. 

VII.  Awards  and  Honors 

Xiaoyan  Yan,  the  graduate  student  supported  during  the  first  year  by  this  project  won  the  outstanding 
graduate  student  award  for  the  NA&ME  Department  for  2008.  Professor  Vlahopoulos,  who  is  the  PI  of 
this  project,  received  the  2011  NA&ME  Departmental  Faculty  award. 

VIII.  Leveraged  Funding 

A  project  was  funded  by  NASA  Langley  for  extending  the  composite  EFEA  developments  for 
rotorcraft  applications,  and  for  validating  the  method  through  comparison  to  test  data. 


Appendix  A  -  Derivation  of  Element  Level  Matrices  for  EFEA  composites 


A  Spectral  Finite  Element  Method  (SFEM)  is  used  for  evaluating  the  group  speed  cg  and  the 
damping  77*  due  to  the  effectiveness  of  SFEM  for  taking  into  account  the  layer-wise  composition 
of  composite  panels  along  with  the  shear  deformation  effects.  The  c*g  and  rj*  values  computed  by 
the  SFEM  are  used  in  the  EFEA  formulation  for  deriving  the  element  level  matrices  of  the 
composite  materials.  The  basic  steps  of  the  SFEM  approach  are  presented  here.  Figure  A.  1 . 
depicts  a  plane  wave  propagating  in  the  positive  x  direction  with  frequency  co  and  wavenumber 
k  in  a  multilayer  composite  plate.  The  x  direction  comprises  a  reference  for  orienting  the  various 
layers.  The  through-thickness  discretization  in  the  z-axis  along  with  the  assumption  of  the 
harmonic  wave  motion  in  the  x  direction  results  in  the  following  form  of  displacement  field, 
ur  =  [u,  v,  w]  at  any  point,  (x,  y,  z )  within  the  plate, 

u(x,z,  t)  =  N(z)ue1<-a,t_fcx-)  0) 

where  N(z)  is  a  matrix  of  shape  functions  and  u  is  a  vector  of  the  nodal  degrees  of  freedom  of 
the  form: 

u=[u1,v1,w1  u2,v2,w2  •••  uNe+1,vNe+1,wNe+1]T  (2) 

where  Ne  is  the  total  number  of  linear  finite  elements  through  the  thickness.  It  is  noted  that  any 
of  the  displacement  components  in  u  may  be  complex  numbers. 


Figure  A.  1 :  A  plane  wave  propagation  in  a  multilayer  composite  panel 
The  time-averaged  total  kinetic  and  potential  energies,  (T)  and  (F)  are  given  by 

=  Ti  (V)  =  \  Jn  oHedn  (3) 

where  H  stands  for  the  Hermitian  transpose.  Substitution  of  appropriate  stress-strain  and  strain- 
displacement  relations  and  the  displacement  field  equation  (1)  into  equation  (3)  results  in: 

(T)  =  uhw2Mu;  (V)  =  u"Ku 

(4) 

where  K  and  M  are,  respectively,  the  stiffness  and  mass  matrices.  The  replacement  of  the  spatial 
derivatives  with  respect  to  x  and  z  by  -  ik  and  dN/dz  can  yield  the  expressions  for  the  stiffness 
and  mass  matrices.  Using  Hamilton’s  principle  results  in: 


(5) 


[K(fc)  -  o>2M]u  =  0 


Since  K  =  K(/c),  for  a  prescribed  circular  frequency  co,  equation  (5)  can  be  solved  for  the 
eigenvalues  co2  and  eigenvectors  u.  The  hysteretic  damping  model  can  be  applied  to  each  layer 
to  yield  the  following  form  for  the  total  time-averaged  energy  loss  associated  with  an  arbitrary 
wave  type, 


NP 
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where  rh  and  (e)l  are,  respectively,  the  structural  loss  factor  and  energy  density  of  the  1th  layer 
of  multi-layered  composites.  Since  (e)  =  (T)  +  (I/)  =  2(1/),  the  damping  loss  factor  associated 
with  each  propagating  wave  can  be  expressed  as 


uHKu 


(7) 


Here,  the  distinction  should  be  made  that  K®  is  the  stiffness  matrix  of  the  1th  layer  and  K  is  the 
assembled  global  stiffness  matrix  of  all  layers.  Equation  (7)  is  derived  for  each  wave  heading 
angle  6h  and  thus  the  angle-average  must  be  considered  for  accounting  for  the  full  range  of  wave 
propagation  directions  when  deriving  the  angle  averaged  damping  loss  factor  as  follows: 


/9  rlWddO, i 
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The  variation  of  the  angle  is  obtained  by  rotating  the  relative  orientation  of  the  x-axis  with 
respect  to  the  composite  material  (Figure  A.  1).  The  circular  frequency  for  the  wave  number  k 
can  be  written  in  terms  of  the  Rayleigh  quotient  as 

uHK(k)u  (9) 

"  =  u"Mu 

Figure  A.2  shows  the  wavenumber  of  an  arbitrary  wave  propagating  in  an  anisotropic  medium  as 
a  function  of  wave  heading,  0;.  The  distribution  and  the  direction  of  energy  flow  in  anisotropic 
media  like  composites  exhibit  directional  dependence,  thus  the  direction  of  the  energy  flow, 
depicted  as  9e  in  Figure  A.2,  is  usually  different  from  the  wave  heading,  0j.  Thus,  an 
appropriate  correction  must  be  introduced  in  the  calculation  of  the  group  speed. 


Figure  A.2:  Wavenumber  as  a  function  of  wave  heading  in  the  k  (wave  number)  plane 

The  group  speed  in  the  direction  of  wave  propagation  is  cg9  =  dw/dk,  and  the  differentiation  of 
equation  (9)  with  respect  to  the  wavenumber,  k,  results  in: 


cg6 


uH  (5K(/c)  /  dk)u 

2muHMu 


Referring  to  Figure  A.2,  the  following  relationship  holds: 

=  _Cg9__ 
9  cos (0£  -  9e) 


(10) 


(11) 


where  the  heading  of  group  speed,  9e  may  be  derived  from  the  geometric  interpretation  of  the 
wavenumber  curve  shown  in  Figure  A.2 


dkx/d0i  (dk(9i)/ddi)cos9i  —  k(9i)smdi  (12) 

dky/ddi  (d/c(0i)/d0j)sin0j  +  /c(0j)cos0j 


Therefore,  the  angle-averaged  group  speed  can  be  evaluated  for  the  full  range  of  wave 
propagation  directions  as  follows: 


J0  Cg(Gj)dGi 

fe  dGt 


(13) 


In  the  EFEA  formulation  for  composites,  c*g  and  r]*  from  Equations  (13)  and  (8)  respectively,  are 
used  when  deriving  the  element  level  EFEA  matrices. 


Appendix  B  -  Derivation  of  Power  Transfer  Coefficients  between  Composite 
Members  and  Derivation  of  Joint  Matrices 

The  derivation  of  power  transfer  coefficients  is  based  on  considering  a  connection  between  semi- 
infinite  members,  prescribing  an  impinging  wave  type  originating  from  one  of  the  members,  and 
evaluating  all  transmitted  or  reflected  wave  types  in  all  members.  Once  available,  the  power 
transfer  coefficients  are  used  for  computing  the  EFEA  joint  matrices  which  connect  elements 
across  a  discontinuity  in  the  EFEA  system  of  equations.  The  formulation  for  computing  power 
transfer  coefficients  between  composite  members  is  presented  here.  It  must  be  noted,  that  this 
formulation  accounts  for  First  order  Shear  Deformation  Theory  (FSDT)  effects,  which  are 
important  when  considering  sandwich  composite  panels.  The  extended  wave  dynamic  stiffness 
matrix  approach  is  used  to  calculate  the  wave  power  transfer  coefficients  of  coupled  composite 
members. 

Figure  B.l.  presents  a  junction  consisting  of  N  semi-infinite  composite  members.  The 
coordinate  systems  and  the  translational  and  rotational  displacements  are  also  presented  in  Figure 
B.  1. 


z,  iv 


(b) 

Figure  B.l:  (a)  A  general  N-plate  junction  and  global  coordinate  system;  (b)  local  coordinate 

system  and  displacements  for  plate  n 


For  either  a  laminated  or  a  sandwich  composite  member,  the  FSDT  allows  writing  the  dynamic 
equations  of  motion  in  the  following  matrix-vector  form: 


Lu  =  0 


(l) 


T 

where  u  =  [u,  v,  w,  cpx,  cpy ]  and  L  is  a  linear  differential  operator  which  can  be  expressed  in 
terms  of  elastic  constants  Atj,  and  and  inertial  properties  Ia(a  =  0,1,2).  Substitution  of  a 
wave  form,  u  =  u e-ikxx+Ay+icot  jnto  equation  (1)  yields 

(A2b2  +  dbx  +  b0)u  =  0  (2) 

d  d  d 

where  b0,  bl5  and  b2can  be  obtained  from  L  by  replacing  — ,  —  and  —  by  (-i/cx),A,  and  co.  This 

quadratic  eigenvalue  problem  can  be  solved  for  five  pairs  of  A’s  and  u’s  for  given  kx  and  co. 

Then,  the  edge  displacements,  dy  and  the  corresponding  elastic  tractions,  ty  of  each  plate  can  be 
expressed  in  terms  of  the  reflected/transmitted  wave  amplitudes  c  =  [cx,  c2,  c3,  c4,  c5]T  and  an 
incident  wave  amplitude  c0  as  follows 

dy  =  [dyl  dy2  "•  dy5]c  +  dyQC() 

ty  =  [tyl  ty2  •"  ty5]c  +  ty0C0  (4) 

where  dyiand  tyi(/=l,2,...,5)  represent  edge  displacements  and  elastic  tractions  due  to  the  ith 
transmitted  or  reflected  wave  and  dy0  and  ty0  are  those  for  incident  waves.  Equations  (3)  and 
(4)  can  be  combined  to  give 

ty  =  Kdy  -  fy0  (5) 

where  K  =  [tyi  ty2  •••  ty5][dyl  dy2  •••  dy5]_1  is  the  wave  dynamic  stiffness  matrix  and 
fyo=  (Kdy0- 1  y0)co  is  the  force  due  to  the  incident  wave.  It  should  be  noted  that  fy0  =  0  in  the 
plates  where  no  incident  wave  exists.  By  implying  subscript  “y”  for  the  remaining  discussion, 
equation  (5)  can  be  rewritten  for  the  nth  plate  as 

tn  =  Kndn  -  fn  (6) 

The  wave  dynamic  matrix  is  assembled  to  yield  an  equation  for  the  computation  of  power 
transmission  coefficients.  Since  dn  and  tn  are  defined  in  the  local  coordinate  system,  they  need 
to  be  transformed  to  those  at  the  common  junction  defined  in  the  global  coordinate  system,  dy 
and  ty  as  follows: 

dn  =  Rndy;  ty  =  Rntn  (7) 

where  Rn  is  a  simple  coordinate  transformation  matrix  consisting  of  cosi pn  and  sini pn.  Since  the 
force  equilibrium  equation  for  the  N  plates  connected  through  a  structural  joint  is  £n=i  R^tn  = 
0,  using  equations  (6)  and  (7)  results  in: 


r N  i 

r N  l 

Y  RnKnRn 

dy  = 

Y  Kfn 

-n=l 

-n=l 

(8) 


where  d j  can  be  evaluated  from  Equation  (8).  Using  equation  (6)  the  wave  amplitudes  are 
evaluated: 


C  —  [dyl  dy2  "•  dy5]  ^(ly  —  dy0CQ) 

Then,  the  power  transfer  coefficients  can  be  readily  calculated  from: 


T;/(m,0)  = 


(9) 

(10) 


where  c;-  is  the  amplitude  of  a  reflected  or  transmitted  wave  type  j  for  the  given  coi,  which  is  the 
amplitude  of  an  incident  wave  type  i  with  frequency  to  and  an  incidence  angle  of  0. 


The  above  computed  transmission  coefficients  are  highly  dependent  on  the  angle  of  incidence,  9 
for  anisotropic  media  (i.e.  composite  panels).  Thus,  the  power  transfer  coefficients  are  obtained 
through  averaging  of  the  power  transmission  coefficients  with  respect  to  the  angle  of  incidence: 


Tjj(,@)cgy(9)  I 

mCge{9 )  / 


vgy 


0 9 ) 


c(9)cge(9 ) 


d9 


(ll) 


The  power  transfer  coefficients  are  used  for  creating  the  EFEA  joint  matrices  which  connect 
elements  across  a  discontinuity.  The  transmission  coefficient  associated  with  each  one  of  the 
generated  waves  is  calculated  as  the  ratio  between  the  power  of  the  transmitted  wave  divided  by 
the  power  of  the  wave  incident  to  the  junction.  The  complete  set  of  transmission  coefficients  for 
the  junction  is  written  in  the  form  T,Jpr  ( co ,  (j>)  where  i,  p  to,  and  ())  represent  respectively,  the  carrier 

plate,  the  wave  type,  the  frequency  and  the  heading  of  the  incident  wave,  and  j ,  r  represent 
respectively  the  carrier  plate  and  the  wave  type  of  the  transmitted  wave.  The  matrix  of  the 
power  transfer  coefficients  is  computed  [r]( .  The  joint  matrices  in  the  EFEA  formulation  define 

the  power  transfer  across  elements  at  the  joints  and  are  derived  from  the  averaged  power  transfer 
coefficients: 

vi = m-vixm+vir'lMjdB  02) 

where  <(); ,  (j^  are  Lagrangian  basis  functions,  B  is  the  boundary  area  between  elements  i  and  j  at 
the  joint,  and  [r]';  is  a  matrix  comprised  by  the  averaged  power  transfer  coefficients  through  the 
connection. 


